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INTRODUCTION 

In the fuel cell power section, air, in excess of the stoichiometric 
mixture, enters the cathode side of the cell, and effluents from the low 
temperature shift converter enter at the anode. The anode input contains 
HgO, H 2 » CO and C0 2 - Jn this analysis, it is assumed that a 
fixed percentage of hydrogen is consumed at the anode, and the H 2 0 being 
formed exits the fuel cell, with the depleted air, through the cathode exit. 
The overall reaction in the fuel cell power section is, 

H 2 + 1/2 0 2 * H 2 0 

Two distinct mathematical models of fuel cells have been developed with 
computer programs for performing the necessary calculations. The first was a 
"lumped parameter" model; the second was a three-dimensional detailed model of 
the stack. 

The simplified lumped model, described in the previous report, is an 
"input-output" model developed for the system trade-off studies (Ref. 1). 

The detailed distributed model is a finite-difference model of the opera- 
tion of the fuel cell which was used to calculate the effects of the cell and 
module design on performance. It calculates the current density distribution 
in the cells as a function of the local reactant compositions, local tempera- 
tures, catalyst utilization factors, etc. Since these are interdependent 
(e.g., the local temperature depends on the local current density), the 
computations are highly iterative and require considerably more computer 
capacity and time than the lumped model. An associated computer program will 
be used to compare an alternative design of cooling scheme in the stack. 


I. LUMPED MODEL AND VOLTAGE-CURRENT CHARACTERISTIC 
1.1 Mass and Energy Balances for Lumped Model 

The lumped model provides a rapid (in terms of computation time) means of 
calculating the fuel cell module output characteristics (voltage, current, and 
heat generation rate) in terms of the inputs from the fuel processing subsystem 
and the gross fuel cell design parameters such as catalyst loading. 


The mass balances of hydrogen, oxygen and water are as follows: 



NX ^2 = NI ^2 - (Imean A)/(nJ) 

(1-1) 


NX.Q 2 = NI 02 - (Imean A)/(2n?) 

(1-2) 


NX H 20 = NIH20 + (Imean A)/(nJ) 

(1-3) 

where NX: 

exit flow rate of hydrogen, oxygen, or steam, g-mole/sec 


NI: 

inlet flow rate of hydrogen, oxygen, or steam, g-mole/sac 


Imean: 

p 

mean current density, A/cm 


A: 

2 

effective area of cell plate, cm 


n: 

number of Faraday equivalents transferred 


7 "- 

Faraday constant 


The energy balance for the fuel cell is 

«(Q + W ) = n o ^ ^ kf)j ”2Z n i ( A 

e' pp J J rF 



(0P) S dT-Znif 298 (Op), dT 
PF J J 298 2 rF Jl ip 

(1-4) 


where the subscripts PF, rF represent the products and reactants in the fuel 
cell, respectively. TfF is the final temperature of the products and TiF is 
the initial temperature of the reactants in the fuel cell. The nj and ni are 


the species flow rates of the products and reactants, respectively. The terms 
Q and W are the rates of heat and the electrical energy generation by the fuel 
cell, respectively. Q is proportional to the specific heat generation Q F 
where: 


Q = N p Xn Yn Qp 
and 0 F = (tt - V) I 


(1-5) 

(1-6) 


where Q: 
Q f : 

V 

Xn: 

Yn: 

I: 

AHr: 


total heat generated, 0/sec 

O 

heat generated per unit area of cell, J/sec car 

number of cells 

width of cell plate, cm 

length of cell plat.'?, cm 

2 

fuel cell current density, A/cm 
heat of reaction, J/g-mole of H 2 


1.2 Voltage-Current Characteristics 

Because of the irreversibility, the voltage V for a working fuel cell is 
the difference between the open circuit voltage and the cell polarization 
terms: 


V = E - n 


(1-7) 


where E: Nernst potential (reversible open circuit E.M.F.) 

n: overpotential or polarization 

The reversible cell potential, E is given by the Nernst equation: 


3 


(1-8) 


7 


RT 


E o - E ( T ) + nr ln 


YH 0 ,r PtY0- 

Cj t 

“iro 

2 


with Pt ; total pressure, atm 

E 0 (T): standard E.M.F. of cell at temperature T, volts 


E 0 (T) = 1.261-0.00025 T, T, K (Ref. 2) 


YH 


2‘ 
Y0 2 : 


mean mole fraction of hydrogen at anode 
mean mole fraction of oxygen at cathode 
YH 2 0: mean mole fraction of water vapor at cathode 

The polarization term n consists of four components, 

n = na + nr + nd + nco 
activation polarization at cathode, volts 
resistance polarization, volts 
diffusion polarization, volts 

activation polarization at anode due to CO poisoning of 
catalyst, volts 


where na: 
nr: 
nd: 
nco: 


and 


with oCo: 
i: 
io: 
SA: 
CL: 
CU: 


a _ Kl In 1 

na “<*oZF ( io) (SA) (CL) (CU ) 
transfer coefficient 
current density, mA/cm 2 
exchange current density of cathode, mA/cm 2 

p 

specific catalyst surface area, cm /g 
catalyst loading on cathode, g/cm 2 
catalyst utilization factor 


(1-9) 


(1-10) 


* 
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The exchange current is a function of the acid concentration, temperature, 
and partial pressure of the oxygen. The acid concentration is a function of 
the water vapor partial pressure which permits correlaion of io as a function 
of Y02, YH20, and T. An empirical fit is 

io « 232.7 (PtY02) 0 * 8 ( PtYH20 ) 0 * 4377 exp (-6652/T) (1-11) 

The resistance polarization is 

nr » ir 

o 

where r: specific cell resistance, ohm-cm . 


The expression of nco was chosen to have strong temperature dependence, 
be directly proportional to Yco, and have a logarithmic dependence on i, iao, 
and catalyst effective area. The resulting expression (Ref. 2) is 


nco = 0.0782PtYco exp 


9 1 9 ° ( j - 4 g 0 ) ( In CL - i 


ao 


( 1 - 12 ) 


where CLa: anode catalyst loading, g/cnr 

r 

iao: anode exchange current, mA/cnf 1 


Diffusion polarization has been neglected here because it is significant 
only at very high current densities. 

In the associated computer code, Subroutine VI, calculates cell voltage as 
a function of the current density or alternatively solves the nonlinear 
equation to evaluate current density as a function of the cell voltage. 


II. CURRENT DENSITY DISTRIBUTION 

In the fuel cell module, the combined modeling of temperature and current 
distribution is an absolute condition for reliable scaling-up of the results 
obtained with small cells, and for predictive models starting from elementary 
porous-electrode representations. 

This subsection describes the calculation of the current density 
distribution over a cell plate on which the air arid fuel flows are at right 
angles. The procedure divides a cell plate into "grids" which are small 
enough so that variations in fuel and oxidant composition and temperature are 
negligible. Then by means of calculation of the boundary conditions for each 
"grid" and iteration, a solution will be obtained that satisfies the input 
specifications (e.g., average current density, fuel and air utilization, and 
reactant flow rates). A diagram of the "grid" is shown in Figure 1. 

The overall method is to first specify a desired average current density i 
for the whole plate and then determine the corresponding voltage V for the 
plate. This voltage will be determined such that it produces unique local 
current densities over the plate whose average value approximates i within a 
specified tolerance. A trial-and-error procedure is used to estimate the 
local current density and overall voltage. The model basically applies the 
same voltage-current equation used in the lumped model {described in Chapter 
1) to each grid section of the cell. 
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N1UH2(1,1) N1UH2(1,2) I N1UII2(l t NY) 



IJ2U02(1,1) N2TJ02(2 f l) 


Oxidant 


Figure 1 Finite Difference Model Definition of Current 
Density Distribution on Cell Plate 





Mathematic Formulation 

Exit flow of hydrogen from gird (ij) 
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NX H 2 (i,j) » NI H 2 (1,j) - (I(U)A)/(n7) (2-1) 

Exit flow of oxygen from gid (i,j) 

NX 0 2 (i ,j) = NI 0 2 (i,j) - (1(1 »j)A)/(2n 7*) (2-2) 

Exit flow of water from grid (ij) 

NX H 2 0(i,j) . NI H 2 0(i,j) + (I(i ,j)A)/( n?) (2-3) 


where NX H 2 , 0 2 , H 2 0(i,j) : hydrogen (oxygen or water) portion flow 

rate at exit of grid (i,j), g-mole/sec. 


NI H 2 , 0 2 , H 2 0(i,j) 


I ( i » J ) 
A 


: hydrogen (oxygen or water) portion 
flow rate at inlet side of grid (i,j), 
g-mc'l«/$ec. 

; current density of grid (i,j), A/cm 2 
: area of grid, cm 2 


The flow charge of executive program (CUPRO) for calculating current density 
distribution is shown in Figure 2. 


CUPRO 
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calculate 
inlet flow ra ;e 
of component 


calculate 
aver, molar 
fractions of 
components 


estimate init: al 
guess of V 
CALL VI^ 



DO I * 1, NX 
DO J * 1,NY 


calculate 

^ VI 1 ri + rnrtl o-w 

V ■■■• ■ ■- 1 

fractions of 
components 


initial guess 
of C.D. 


calculate 
outlet flow 
rate 


calculate 
outlet molar 
fractions of 
components 



correct 
initial 
guess of 
C.D. 


estimate 
C.D. of (I,J 
CALL VI^ 


VA 

^/test^N. 

-&< 

initial gue^s 


\C.D^/ 


YES 


nnurnTOrra 

VW*»A 0m 


calculate 
avg. C.D. 
for plate 


^design ref 
^N^vg. CjJU 


WRITE 

the 

results 


correct 
initial 
guess of 
V 


RETURN 


(1) : M=1 yields voltage 

(2) : M=2 yields C.D. 


Figure 2 Flow Chart of CUPRO 


















III. THERMAL ANALYSIS AND TEMPERATURE DISTRIBUTION 

The electrical energy proauction in phosphoric acid fuel cells is 
accompanied by approximately equal amounts of heat energy generation. Removal 
of this heat can be accomplished by a suitable flow of input gases or by using 
separate cooling plates. 

The work reported in this section is directed towards estimating the 
steady state temperature profiles in practical phosphoric acid fuel cell 
stacks. The fuel cell stack considered in this section is composed of cell 
plates on which the air (oxygen) and fuel (hydrogen) flows are at right 
angles. A cooling plate is placed between individual groups of cells at a 
regular interval. Symmetry in the stacking direction occurs at the middle of 
a cooling plate and midway between cooling plates. 

3.1 Previous Work 
«' 11 — — ■ - ■■■■■'■■ 

Estimation of the temperature profiles in an operating cell is important 
for the estimation of the power density distribution, thermal stability, and 
cooling requirements. Only a limited amount of information on this subject 
has been reported in the past. Baker and coworkers recognized this need and 
have performed a comprehensive study of steady state heat transfer in 
electrochemical systems (Refs. 3, 4, 5). They studied various cases involving 
one dimensional analysis of a single adiabatic fuel cell and a three 
dimensional analysis of a multicell stack. 

A single fuel cell with no lateral heat transfer and no conduction of heat 
through the cell in the direction perpendicul ar to the gas flow was considered 
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(Ref. 4). Heat transfer by conduction in the direction of the gas flow was 
considered negligible in comparison to the heat transfer by convection, and 
analytical expressions for the electrolyte, fuel, and air temperature profiles 
were derived. 

For the three dimensional analysis of the stack, it was assumed that all 
of the walls except for the wall from which the air enters were maintained at 
a constant temperature. The rate of heat generation per unit volume of the 
stack was assumed constant. An analytical solution for the temperature 
profile was developed, assuming that the electrolyte and gas temperatures were 
not very different. 

Another paper (Ref. 5) considered various limiting and special cases to 
determine the maximum temperature of a stack. Two dimensional heat transfer 
analysis was carried out in the case of a thick stack where heat transfer in 
the direction of stacking was neglected. In the case of thin stacks, three 
dimensional heat transfer was considered with each wall at a different 
temperature. Infinite series solutions were developed for both thick and thin 
stacks. The authors estimated the maximum stack temperature for the constant 
wall temperature case. An approximate formula to predict the effect of 
conductivities, size, and current density on the maximum stack temperature was 
developed. A generalized analysis, which can incorporate the effect of finite 
resistance to heat transfer at the wall, the effect of cold or hot feeds, or 
nonuniform heat generation, was also carried out using the method of Green's 
function. 
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3.2 Temperature Distribution 

The temperature distribution for the module was developed from the 
temperature distributions within representative slices or strips within a set 
of cel v and cooling plate cells. The analysis includes conduction within 
bipolar plates, conduction between plates, the separate cooling effects of the 
process air and the coolant (basically air is considered as the coolant), and 
the temperature change of air flows along their respective channels. The 
distribution of the heat generation is determined from the current density 
distribution. 

The model assumes that (1) the temperature gradients in the direction of 
the fuel flow are small. This assumption is justified since the major 
temperature gradients are in the air flow direction and since the heat 
capacity of the fuel stream is only a few percent of the heat capacity of the 
air stream; (2) the edge of the cell is operating adiabatically; (3) a half 
set of cell plates between cooling plates is analyzed, which includes one half 
cooling plate and two and a half ceil plates. Thus, because of the symmetry, 
all of the stack behaves similarly. The geometry of a representative slice 
(Lx x Ly x Lz) through the stack is shown in Figure 3. 

Mathematical Formulation 

The material balances of the fuel and the oxident have been presented in 
Chapter 2. There are four energy balance equations for the cell plate, 
cooling plate, process air, and coolant. 
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ORIGINAL PAGii & 
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cel 1 on process air side in air flow direction 

dl I 


t Kv/® 2t + K* 3T 

1 Ky TJ Kx ax 


- Kx — I - — E- d JR ( v*-Vl I 

x+t NX ax i x pp ay " ' } . 


(3-1) 


cooling plate in coolant direction 


f Ky ^ + 2Kx ?T 

a y 


dx 


Cc m 


c d Tc 


x + t'/2 Pc ay 


= 0 


process air side 


coolant side 


±lR = h P , S ip . (T _ T j 
d y mp Cp u 1 p } 


d Tc _ he Sc 
d y ” me Cc 


(T-Tc) 


Boundary conditions 
x = 0 aT/ax = 0 


y - 0 

x = Lx 
y = Ly 
y = 0 

y - o 

where m 
C 

Ky 

Kx 

t 

xl 

x2 


syrranetric condition 
adiabatic assumption 
symmetric condition 
adiabatic assumption 


(3-2) 


(3-3) 


(3-4) 


9T/ay = 0 
aT/^x = 0 

aT/ay = o 

Tp = Tp, inlet 
Tc = Tc, inlet 

= mass flow rate, Kg/hr-channel 
= heat capacity, J/Kg-K 

= effective thermal conductivity of cell in flow direction, 
J/hr-m-K 

= effective thermal conductivity of cell on stacking direction, 
j/hr-m-K 

= thickness of cell including fuel and air channel, m 
» effective conduction distance from plate to upper cell plate, m 
= effective conduction distance from plate to lower cell plate, m 
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MU , a 


ORIGINAL L£] 

OF POOR QUALITY 

P « pitch of channel , m 

xl ' = effective conduction distance from cooling plate to upper cell 

plate, m 

Lx, Ly = height and length of one slice, respectively, m 
V* = aH/ZF, V 

t* = thickness of cooling plate, m 

h = heat transfer coefficient, J/hr-m^K 
S = perimeter of the channel, m 
Subscription 

p .- process air 

c = cooling air 

These simultaneous ordinary differential equations and corresponding boundary 
conditions were solved by the finite-difference method. The final difference 
equations are in next subsection. 

Finite-Difference Model 

The energy balance on an internal element j (2<j<N~l) for bipolar plate i 
(2<i<Nl) can now be written as (sea Figure 3) 

- (^)Ti,j-l+(2 i^4 + JT + vf) Ti, J 
aY^ aY^ * a 

" ( ^2^ Ti,y+1- (|j) Ti-l,j-(||) Ti+l,j (3-5) 

+ ^^aY^ ( Tpi ’J " Tpisj “^ = ^*~V) 11 

The energy balance on an internal element j (2<_j<N-l) of the cooling plate i=l 
can be written as 
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O' ■' 
* 





I 

i 


/ 

,? 

J 



orssal PX3 E2T 

OF POOR QUAUifi? 

- (^-r-) T1, 0-1 + (2 + 2 |p-) Tl, j 


aV 


aY 


(**-§") Tl, j+1 - (2 §•) T2, j 

&r w 


+ (TCJ-TCj-1) . 0 


(3-6) 


The energy balance on interior element j (2<_j<N-l) of the symmetric plate i=Nl 
is 


- <Jr> t ni,o-i + ‘ 2 Jfr + tj) T w,j 


- t 


aY 


2 ' Nl,j+1 v X2' N 1— 1 , j 


/ 6,Kx x T 

v Y 9 < U 


(3-7) 


+ O (T » N i,r Tp Ni,j-i) ■ < V *-V ‘nuj 

2 for odd values of NK 
1 for even values of NK 

the number of cell plates between two adjacent cooling plates 
1 + NK/2 for even NK 
1 + (NK+U/2 for odd NK 


where 6 = 

5 = 

NK: 

Nl: 


The energy balance on elements j=l are obtained as above, except for: the 

values of Ti,0 are replaced by Ti,l: the values of Tpi,0 are replaced by TPO, 
which is the inlet process air temperature; and TCo is replaced by TCO, the 

V 

inlet cooling air temperature. The energy balances on elements j=N are 
obtained from the above with Ti,j+1 replaced by Ti,N. 
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For the process air flow, one can set up NlxN equations of the form 

TPi,j - TPi , j-1 + (Ti,j - TPi ,j~l) (l-e^ pi »j) (3-8) 

where 

hi , SP 

^pijj " Mp Cp aY ' 

For the cooling air flow, one obtains N equations of the form 

TCj = TCj-1 + (Tl,j - TCj-lHl-e"^) (3-10) 

where 

. h iSc 

ij ■ « ( 3 - U > 

Thus, the total number of temperature equations matches the number of unknown 
temperatures and the set can be solved using the Gaussian elimination method 
with calculated or input values of cell voltages, current densities, mass 
flow, heat generation and heat transfer coefficients. Each resulting tempera- 
ture distribution is used to recalculate the current density distribution 
until covergence is reached. The relationship between voltage and current and 
the calculation of heat generation have been presented in Chapter 1. 

Heat Transfer Coefficients 

An empirical equation (Ref. 6) for the Nusselt number for fully developed 
laminar flow in a rectangular channel is: 

Nuf = 3.61 + 4.63 (1 ~oC) 3 - 2 (3-12) 

where c>^= a/b; a is the smaller side of rectangular channel and 
b is the larger side of the channel. 
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Near the inlet of a channel, the heat transfer coefficient is larger than the 
fully developed value due to development of the laminar boundary layer. If R 
is the ratio of the average Nusselt number for the region 0 to x to the fully 
developed Nusselt number, then (Ref. 7) 


R , . 0.0183 Gz 

1+0.04 Gz 2/3 

where GZ: Graetz number = Re Pr (D^ x ) 

Re: Reynolds number based on 

Pr: Prandtl number of gas 

D^: Hydraulic diameter, m 


(3-13) 




For turbulent flow, the average Nusselt number over the region 0 to x is 
described as (Ref. 8) 

Nu t = 0.116 [Re 2/3 -125] Pr 1/3 [1 + (D/x) 2/3 ] (3-14) 

i 

The flow chart of the executive program (MAIN program) for calculating the 
temperature distribution in the stack is shown in Figure 4. 
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IV. COMPUTER CODE 
4.1 Program Description 


The computer code contains one executive program (MAIN program) and eleven 
subroutines. The mathematical model and algorithm used in MAIN program was 
shown in Chapter 3. Table 1 lists the nomenclature of the program. 

All of the subroutines are listed in Table 2 associated with their 
specified functions. Among these, Subroutines VI and CUPRO have been 
described in Chapters 1 and 2, respectively. Subroutine DRAWE, which execute 
the contour drawing package, will not be used except running the program on 
IBM 370 of NASA Lewis lisearch Center. The rest of listed subroutines are 
used to estimate the properties of the process fluids or for I/O usage. 

The whole program listing is shown in the end of this manual. 

4.2 Program Operation 

The program input only consists of a set of NAMELIST data in a specified 
order. The first NAMELIST set is called DIMEN and contains the dimensions of 
cell and cooling plates, number of cell plates between two cooling plates, 
number of cell plates between two cooling plates, number of air and fuel 
channels and utilization, pressure, number of finite difference sections, and 
input temperature on anode and cathode sides. The order of input data inside 
one NAMELIST need not be fixed. 
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TABLE 1 


PROGRAM NOMENCLATURE 


AL: 

AL1(I) : 


ALFA: 

AMUC: 

AMWA: 

AMWC: 

CL: 

CU: 

CM(I): 
CMC ( I ) : 
CPC: 


DNSA(I): 
DNSC(l) : 
DX: 


FCONST : 
G(I,J): 
GZ(I): 
GZA: 


H2: 

HC(I) : 
h'H: 

PPRO(I ,J,K) : 
KX: 

KY: 

HA: 

MAC ( I } : 

NC: 

NCA: 


NCC: 

NK: 

NP: 

NX: 

NY: 

02 : 

00 : 

PC: 


PHl(I.J) : 
PH2(I): 


POP: 

PP: 

PR: 

QW(I,J): 


aspect ratios of process air 

aspect ratios of cooling channel in different sections (treed 
form) 

transfer coefficient 
viscosity; Ib/hr-ft 

molecular weight of process air; Ib/lb-mole 
molecular weight of cooling air; lb/ lb-mole 
catalyst loading; mg/cm 2 
catalyst utilization 

mole fraction of component I in cooling air 

mole fraction of component I in process air 

heat capacity; Btu/ lb-mo le-R 

moles of component I in process air; lb-mole 

moles of component I in cooling air; lb-mole 

length of x— division; ft 

Faraday constant; 96500 coul ./g-equivalent 

coefficient of simultaneous linear equations 

Graetz number of different sections in cooling channels 

Graetz number in process air 

heat transfer coefficient of plate I x-di vision 0 y-division K 
of process air; Btu/ft 2 -hr-R 
required hydrogen; g-mole/hr~stack 

heat transfer coefficient of division I in cooling channel 
required hydrogen; g-mole/sec-plate 

current density of plate I x-di v i s i on J y-division K; A/cm 2 

effective thermal conductivity in stacking direction; Btu/hr-ft-R 

effective thermal conductivity in flow direction; Btu/hr-ft-R 

mass flow rate of process air; Ib/hr-chanrsel 

mass flow rate of cooling air in section I; Ib/hr-channel 

number of stoich air in cooling channel 

number of process air channels 

number of cooling channels 

number of plates between cooling plate 

number of plates in a stack 

number of divisions in x direction 

number of divisions in y direction 

required oxygen; g-moles/hr-stack 

required oxygen; g-moles/’ sec-plate 

pitch of cooling channel; ft 

dimensionless group of plate I division J in process air 
dimensionless group of division I in cooling air 
inlet gas pressure; atm 
pitch of process air; ft 
Prandtl number of gas 

heat generation rate of division J plate I; Btu/hr 
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TABLE 1 (cont'd) 
PROGRAM NOMENCLATURE 


R(I): 

RA: 

RE( I) : 

SA: 

SRO: 

T: 

Tl: 

TAIN: 

TKA: 

TCI N : 

TKC: 

TDNSC: 

TDH2 ( I ,J) : 

TDH20 ( I , J ) : 

TD02 ( I ,J) : 

TD1(I,J): 

TD2(I,0): 

TFA: 

TFC: 

TAK: 

TCK: 

TKAA: 

TKCC: 

TKF: 

TRR(I): 

TUN: 

UTA: 

UTH: 

WA: 

WAD: 

WAW: 

WC: 

WCD: 

WCW: 

WE: 

WFD: 

WFW: 

WP: 

X( I) : 

XAMP: 

XDNSCO: 


ratio of average Nusselt number for region 0 to x to the fully 
developed Nusselt number of division I in cooling channel 
ratio of average Nusselt number for region 0 to x to the fully 
developed Nusselt number of division I in process channel 
Reynolds number of division I in cooling channel 
catalyst surface area; cm 2 /mg 
cell resistance at 450 K; Qhm-cm 2 
thickness of cell including process channels; ft 
thickness of cooling plate; ft 
inlet temperature of process air; R 
inlet temperature of process air; K 
inlet temperature of cooling air; R 
inlet temperature of cooling air; K 
total moles in cooling channel; g-mole/hr-division 
flow rate of hydrogen in fuel channel at division J plate i; 
g-mole/sec 

flow rate of water in process air channel at division J plate I; 
g-mole/sec 

flow rate of oxygen in process air channel at division J plate I; 
g-mole/sec 

total flow rate in fuel channel; g-mole/sec 

total flow rate in process air channel; 9 -mole/sec 

inlet temperature of process air; F 

inlet temperature of cooling air; F 

thermal conductivity of process air; Btu/hr-ft-R 

thermal conductivity of cooling air; Btu/hr-ft-R 

average temperature of process air; K 

average temperture of cooling air; K 

inlet temperature of fuel; K 

average operating temperature of plate I; R 

Nusselt number 

utilization of air 

utilization of fuel 

hydraulic diameter of process air channel; ft 

depth of process air channel; ft 

width of process air channel; ft 

hydraulic diameter of cooling channel; ft 

depth of cooling channel; ft 

width of cooling channel; ft 

thickness of cell; ft 

depth of fuel channel; ft 

width of fuel channel; ft 

thickness between two cooling plate; ft 

solution of simultaneous equations 

amp/plate 

designed current density; amp/cm 2 
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TABLE 1 (cont'd) 
PROGRAM NOMENCLATURE 


XN: 

XXOO( J,K): 
Yl: 


Y2: 


Y1CH 4 ; 

Y1CO: 

YlCOp : 

Y1H 2 : 

YlHpO : 

Y2H 2 0: 

Y2N 2 : 

Y20 2 ; 

YN: 

Z: 


length of cell in x-direction; ft 
same as PPRO(I,J,K) in each plate; A/cm 2 
effective conduction distance from cell plate to cooling plate; 
ft 

effective conduction distance from cell plate to another cell 
plate; ft 

mole fraction of CH 4 in fuel 

mole fraction of CO in fuel 

mole fraction of C0 2 in fuel 

mole fraction of H 2 in fuel 

mole fraction of H 2 0 in fuel 

mole fraction of H 2 0 in air 

mole fraction of N 2 in air 

mole fraction of 0 2 in air 

length of cell in ^-direction; ft 

number of Faraday equivalents transferred 


TABLE 2 



DEFINITIONS OF SUBROUTINES 

Subroutines 

DESCRIPTION 

DATAIN 

X. input data reading 

2. changing units 

3. calculation of the constants used in MAIN program 

DATACA 

calculations of the properties and coefficients for 
cooling air 

VI 

calculation of the relationship between voltage and 
current density for specified fuel cell plate 

CUPRO 

estimation of the steady state current density 
distribution on the cell plate 

GAUSS 

Gauss-Seidle iteration used to solve simultaneous 
linear equations 

CMASS 

calculation of the mass fraction of gas stream 

CMOLE 

calculation of the mole ^action of gas stream 

HTCP 

estimation of the heat capacity of specified gas 
mixture 

THC 

estimation of the thermal conductivity of specified 
gas mixture 

VIS 

estimation of the viscosity of specified gas mixture 

ORAWE 

execution of the contour drawing package 
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The second set (ERR) only contains the convergence criterion for program 
trial-and-error procedure. The third NAMELIST set (CZ) specifies the kinetic 
data of the catalyst used in anode and cathode sides. 

PIGA carries the information of coolant flow rate, the dimension of 
cooling channels, and the thermal conductivites along flow direction and stack 
direction. 

The last NAMELIST set contains the inlet compositions of both anode and 
cathode sides. 

All of the input variables are listed in Table 3, along with thier units 
and numerical values in the sample run, which will be discussed in the next 
chapter. 
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TABLE 3 

INPUT DATA FOR 3-D C.P. AND TEMPERATURE DISTRIBUTIONS (STEADY STATE) 


NAMELIST 

LIST 

VARIABLE 

NAME 

SAMPLE 

VALUE 

UNIT 

DIMEN 

XN 

17 

in 

DIMEN 

YN 

12 

in 

DIMEN 

DNSCO 

0.325 

A/ cm^ 

DIMEN 

UTA 

0.5 


DIMEN 

UTH 

0.75 


DIMEN 

POPC 

3.4 

atm 

DIMEN 

POP 

3.4 

atm 

DIMEN 

TKA 

443 

K 

DIMEN 

WFD 

0.00333 

ft 

DIMEN 

WFW 

0.01 

ft 

DIMEN 

NCC 

30 


DIMEN 

WE 

0.00333 

ft 

DIMEN 

TKF 

450 

K 

DIMEN 

T 

0.0108 

ft 

DIMEN 

NK 

5 


DIMEN 

WAD 

0.00333 

ft 

DIMEN 

WAW 

0.01 

ft 

DIMEN 

NP 

23 


DIMEN 

NCA 

80 


DIMEN 

NF 

55 


DIMEN 

T1 

0.02917 

ft 

DIMEN 

NX 

12 


DIMEN 

NY 

12 


DIMEN 

TINGS 

191 

C 

ERR 

ER 

0.01 


CZ 

CLCA 

0.52 

mg/cm 2 

CZ 

CLAN 

0.34 

mg/cnr 

CZ 

CU 

0.15 

CZ 

SA 

500 

cm 2 /mg 

CZ 

SRO 

0.44 

-cm 2 

CZ 

ALFA 

0.5 


CZ 

DKC 

240000 

A/atm 

CZ 

R 


J/ ( g-mol ) (K) 

FALA 

Z 


g-equivalent 

FALA 

DIGA 

FCONST 

NC 


C~g/equivalent 

DIGA 

KX 


Btu/(ft-h-R) 

DIGA 

KY 


Btu/ (ft-h-R) 


' definition 

length of cell plate in x-direction 
length of cell plate in y-direction 
designed current density 
utilization of U 2 in stack 
utilization of H 2 in stack 
pressure of cooling air 
operating pressure in stack 
inlet temperature of process air 
depth of fuel channel 
width >f fuel channel 
number of cooling channels 
thickness of cell (electrode and matrix) 
inlet temperature of fuel 
thickness of cell plate 
number of plates between two cooling 
plates 

depth of process air channel 
width of process air channel 
number of cell plates 
number of process air channels 
number of fuel channels 
thickness of cooling plate 
finite difference number in x-direction 
finite difference number in y-direction 
initial guess of plate temperature 
criterion for convergence 
catalyst loading on cathode side 
catalyst loading on anode side 
utilization of catalyst 
surface area of catalyst 
cell resistance at 450 K 
transfer coefficient 
constant to calcuate limiting current 
density 
gas constant 

number of Faraday equivalents transferred 
Faraday constant 

ratio of cooling air to air consumed in 
stack 

effective thermal conductivity in stack- 
ing direction 

effective thermal conductivity in flow 
direction 


TABLE 3 (cont'd) 


INPUT DATA FOR 3-D C.P. AND TEMPERATURE DISTRIBUTIONS (STEADY STATE) 
NAMELIST VARIABLE SAMPLE 

LIST NAME VALUE UNIT DEFINITION 


DI 6 A 

TKC 


K 

DIGA 

wcw 

0.22 

ft 

DIGA 

WCD 

0.22 

ft 

FUEL 

YlHo 

0.76 


FUEL 

yic6 2 

0.24 


FUEL 

Y1C0 

0 


FUEL 

FUEL 

YICH 4 

YlHoO 

0 

0 


FUEL 

YlNo 

0 


FUEL 

Y20? 

0.208 


FUEL 

Y2Np 

0.782 


FUEL 

Y2H ? 0 

0.01 

lbm/ft 3 

HEATC 

RHOP 

163 

HEATC 

RHOC 

135 

lbm/ft 3 

HEATC 

CCP 

0.25 

Btu/(lbm-R) 

HEATC 

CCC 

0.201 

Btu/{ Ibm-R) 


inlet cooling air temperature 
width of cooling channel 
depth of cooling channel 
mole fraction of H 2 in anode inlet 
mole fraction of C0 2 in anode inlet 
mole fraction of CO in anode inlet 
mole fraction of CHa in anode inlet 
mole fraction of H 2 0 in anode inlet 
mole fraction of N 2 in anode inlet 
mole fraction of 0 2 in cathode inlet 
mole fraction of N 2 in catlioue inlet 
mole fraction of H 2 0 in cathode inlet 
density of cell plate 
density of cooling plate 
heat capacity of cell plate 
heat capacity of cooling plate 


VMA* J v * • • £* 


- jr 


V. SAMPLE PROBLEM 
5.1 Sample Problem 

The distribution of temperature and the accompanied current density 
profiles in the fuel cell stack with 17"xl2" cell plate have been determined 
from the developed computer program. These distributions are shown in numbers 
at each corresponding grid. It is noted that the set of fuel cell stack 
considered is the symmetric part of cell plates between two cooling plates 
(Figure 3). The associated operating voltage of each considered cell plates 
is also shown in numbers. 

The input data, which is discussed in the previous chapter, is displayed 
in Figure 5. Figure 6 contains the output generated by the sample data input, 
where the input data is reprinted first. Next, the operating voltage, the 
current density of each grid, and the temperature of each grid on the cell 
plate numbered from outmost plate to central plate are printed. The last 
piece of information printed is the average operating temperature, the 
operating pressure, and the DC outlet of the specified stack. 

If the program was run on IBM 370 in NASA Lewis Research Center, the 
subroutine DRAWE can be called to draw the contours of different temperature 
levels. Figure 7 shows one of these drawings. 

The CPU time depends quite on the trial-and-error procedure. The initial 
temperature guesses, the criteria of convergence, and the number of finite 
difference sections will determine the computation time. Usually, the CPU 
time to run this code on IBM 370 is about 1 minute. 
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The distribution of temperature and the accompanied current density 
profiles in the fuel cell stack with 17"xl2" cell plate have been determined 
from the developed computer program. These distributions are shown in numbers 
at each corresponding grid. It is noted that the set of fuel cell stack 
considered is the symmetric part of cell plates between two cooling plates 
(Figure 3). The associated operating voltage of each considered cell plates 
is also shown in numbers. 

The input data, which is discussed in the previous chapter, is displayed 
in Figure 5. Figure 6 contains the output generated by the sample data input, 
where the input data is reprinted first. Next, the operating voltage, the 
current density of each grid, and the temperature of each grid on the cell 
plate numbered from outmost plate to central plate are printed. The last 
piece of information printed is the average operating temperature, the 
operating pressure, and the DC outlet of the specified stack. 

If the program was run on IBM 370 in NASA Lewis Research Center, the 
subroutine DRAWE can be called to draw the contours of different temperature 
levels. Figure 7 shows one of these drawings. 

The CPU time depends quite on the trial-and-error procedure. The initial 
temperature guesses, the criteria of convergence, and the number of finite 
difference sections will determine the computation time. Usually, the CPU 
time to run this code on IBM 370 is about 1 minute. 
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Figure 5 continued 
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Figure 6 continued 
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Figure 6 continued 
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Figure 6 continued 
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5.2 Further Developments 


Parametric Sensitivity and Cooling Scheme 

The plate temperature is a function of the current density, the concentra- 
tions of hydrogen and oxygen, and the cooling effectiveness. In order to 
achieve the optimum design with respect to the temperature distribution, more 
studies of the parameters involved and the cooling scheme are necessary. The 
computer model discussed in the previous chapters is used to examine and 
compare these design parameter®. 

The examined parameters include dimension and size of cell plate, thermal 
conductivities in stack and flow directions, average current density, coolant 
flow rate and inlet temperature of process air. 

There are three configurations of cooling channels considered, whose 
nomenclature and definitions are as follows: 

1. Straight: the dimensions of cooling channel are fixed. 

2. Branch: the cooling channel is branched along the coolant flow 

direction, one example is in figure 8. 

3. Varying Width: the width of cooling channel is different along 

the fuel flow direction. 

After the cooling stream has become fully developed the heat transfer 
coefficient drop dramatically. The "branch" configuration was designed to 
prevent the formation of fully developed flow and to increase the flow rate 
(as the total crossectional area is decreased). The "varying width" 
configuration will put more "jolant on the larger heat generation side, but 
the heat transfer coefficient does not change. 
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Figure 8 One "Branc! 


More detailed descriptions and results were shown in References 9, 10, 
and 11. 

Transient State 

In addition, load change is an important and frequent operation in the 
powerplant. Since the PAFC system can be subjected to sudden load changes and 
load ramping, an understanding of the effects of these transient conditions on 
the PAFC system's performance is essential for the optimal design and control 
of the system. The transient change of temperature distribution in the load 
ramping period was simulated by studying the dynamics of the fuel cell stack. 
The results were shown in References 10 and 11. 
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